
logistic <- function (x) exp(x)/(1 + exp(x))
std.err  <- function (x) sqrt(diag(vcov(x)))
rmse     <- function (x, par) sqrt(mean((x - par)^2))
cover    <- function (x, se, par, alpha = 0.05) {
				z <- qnorm(1 - alpha/2)
				mean(ifelse(par >= x - z * se & par <= x + z * se, 1, 0))
			}

# list experiment
generateComparisonData <- function (n, corrs, gammas, deltas, n.control.items, p.t, n.groups) {
  
  data <- mvrnorm(n = n, mu = rep(0, 3), Sigma = matrix(corrs, nr = 3, nc = 3, byrow = TRUE))
  X    <- cbind(1, data[, 1:2])
  control.items <- n.control.items * logistic (X %*% gammas)
  control.items <- round(control.items)

  z <- logistic(X %*% deltas)
  z <- rbinom(n, 1, z)

  t <- rbinom(n = n, size = 1, prob = p.t)

  y <- ifelse(t == 0, control.items, control.items + z)
  
  cutpoints <- quantile(data[, 3], probs = seq(0, 1, length.out = n.groups + 1))
  group <- cut(data[, 3], breaks = cutpoints, labels = FALSE, include.lowest = TRUE)

  data.frame(group, z, y, x1 = data[, 1], x2 = data[, 2], t)

 }


generateComparison <- function (n.sample, data, n.control.items, 
	h, grouping = c("overall", "subgroup"), weighting = "efficient") {
    
	sample.data <- data[sample(nrow(data), n.sample, replace = FALSE), ]
    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
      else group <- sample.data$group
    if (is.null(names(h))) names(h) <- sort(unique(group))

	regular.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
		J = n.control.items, method = "nls"))
	augment.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
		J = n.control.items, method = "nls", 
		h = h, group = group, matrixMethod = weighting))
   
    while(class(regular.fit) == "try-error" | class(augment.fit) == "try-error") {
		sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
	    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
	      else group <- sample.data$group
	    if (is.null(names(h))) names(h) <- sort(unique(group))

		regular.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
			J = n.control.items, method = "nls"))
		augment.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
			J = n.control.items, method = "nls", 
			h = h, group = group, matrixMethod = weighting))
    }
   
    rbind(c(coef(regular.fit), std.err(regular.fit), 0), 
    	  c(coef(augment.fit), std.err(augment.fit), 1))

}


generateMisspecifiedData <- function (n, corrs, gammas, deltas, sq.coef, n.control.items, p.t, n.groups) {
  
    data <- mvrnorm(n = n, mu = rep(0, 3), Sigma = matrix(corrs, nr = 3, nc = 3, byrow = TRUE))
    data[, 1] <- ifelse(data[, 1] > 0, 1, 0)
    X    <- cbind(1, data[, 1:2])

    control.items <- n.control.items * logistic (X %*% gammas)
    control.items <- round(control.items)

    z <- logistic(X %*% deltas + sq.coef * data[, 2]^2)
    z <- rbinom(n, 1, z)

    t <- rbinom(n = n, size = 1, prob = p.t)

    y <- ifelse(t == 0, control.items, control.items + z)
    
    cutpoints <- quantile(data[, 3], probs = seq(0, 1, length.out = n.groups + 1))
    group <- cut(data[, 3], breaks = cutpoints, labels = FALSE, include.lowest = TRUE)

    data.frame(group, z, y, x1 = data[, 1], x2 = data[, 2], t)

 }


generateMisspecifiedComparison <- function (n.sample, data, n.control.items, 
	h, grouping = c("overall", "subgroup"), weighting = "efficient") { 

	sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
      else group <- sample.data$group
    if (is.null(names(h))) names(h) <- sort(unique(group))

	regular.fit <- ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
		J = n.control.items, method = "nls")
	augment.fit <- ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
		J = n.control.items, method = "nls", 
		h = h, group = group, matrixMethod = weighting)
   
    while(class(regular.fit) == "try-error" | class(augment.fit) == "try-error") {
		sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
	    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
	      else group <- sample.data$group
	    if (is.null(names(h))) names(h) <- sort(unique(group))

		regular.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
			J = n.control.items, method = "ml"))
		augment.fit <- try(ictreg(y ~ x1 + x2, data = sample.data, treat = "t", 
			J = n.control.items, method = "nls", start = coef(regular.fit),
			h = h, group = group, matrixMethod = weighting))
    }
    
    sample.X <- cbind(1, sample.data$x1, sample.data$x2) 
    sample.0 <- sample.X[which(sample.data$x1 == 0), ]
    sample.1 <- sample.X[which(sample.data$x1 == 1), ]

    result <- rbind(
    	            cbind(mean(logistic(sample.0 %*% coef(regular.fit)[1:3])), 
    	            	  mean(logistic(sample.1 %*% coef(regular.fit)[1:3])), 
    	            	  0), 
    	            cbind(mean(logistic(sample.0 %*% coef(augment.fit)[1:3])), 
    	            	  mean(logistic(sample.1 %*% coef(augment.fit)[1:3])), 
    	            	  1)
    	            )
    colnames(result) <- c("x1.0", "x1.1", "augment")

    result

 }	


# randomized response experiment
generateComparisonDataRr <- function (n, corrs, deltas, p.t, n.groups) {
  
  data <- mvrnorm(n = n, mu = rep(0, 3), Sigma = matrix(corrs, nr = 3, nc = 3, byrow = TRUE))
  X    <- cbind(1, data[, 1:2])
  
  z <- logistic(X %*% deltas)
  z <- rbinom(n, 1, z)

  t <- rbinom(n = n, size = 1, prob = p.t)

  y <- ifelse(t == 1, 1, z)
  
  cutpoints <- quantile(data[, 3], probs = seq(0, 1, length.out = n.groups + 1))
  group <- cut(data[, 3], breaks = cutpoints, labels = FALSE, include.lowest = TRUE)

  data.frame(group, z, y, x1 = data[, 1], x2 = data[, 2], t)

 }


generateComparisonRr <- function (n.sample, data, h, grouping = c("overall", "subgroup"), weighting = "efficient") {
    
  sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
  if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
	else group <- sample.data$group
  if (is.null(names(h))) names(h) <- sort(unique(group))

  regular.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known"))
  augment.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known", 
  	h = h, group = group, matrixMethod = weighting))
   
  while(class(regular.fit) == "try-error" | class(augment.fit) == "try-error") {
  sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
      else group <- sample.data$group
    if (is.null(names(h))) names(h) <- sort(unique(group))

  regular.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known"))
  augment.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known", 
  	h = h, group = group, matrixMethod = weighting))
  }
   
  rbind(c(coef(regular.fit), std.err(regular.fit), 0), 
        c(coef(augment.fit), std.err(augment.fit), 1))

}


generateMisspecifiedDataRr <- function (n, corrs, deltas, sq.coef, p.t, n.groups) {
  
    data <- mvrnorm(n = n, mu = rep(0, 3), Sigma = matrix(corrs, nr = 3, nc = 3, byrow = TRUE))
    data[, 1] <- ifelse(data[, 1] > 0, 1, 0)
    X    <- cbind(1, data[, 1:2])
    
    z <- logistic(X %*% deltas + sq.coef * data[, 2]^2)
    z <- rbinom(n, 1, z)

    t <- rbinom(n = n, size = 1, prob = p.t)

    y <- ifelse(t == 1, 1, z)
    
    cutpoints <- quantile(data[, 3], probs = seq(0, 1, length.out = n.groups + 1))
    group <- cut(data[, 3], breaks = cutpoints, labels = FALSE, include.lowest = TRUE)

    data.frame(group, z, y, x1 = data[, 1], x2 = data[, 2], t)

 }


generateMisspecifiedComparisonRr <- function (n.sample, data, 
  h, grouping = c("overall", "subgroup"), weighting = "efficient") { 

  sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
    if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
      else group <- sample.data$group
    if (is.null(names(h))) names(h) <- sort(unique(group))

  regular.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known"))
  augment.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known", 
    h = h, group = group, matrixMethod = weighting))
   
    while(class(regular.fit) == "try-error" | class(augment.fit) == "try-error") {
    sample.data <- data[sample(nrow(data), n.sample, replace = TRUE), ]
      if (grouping == "overall") group <- rep("1", nrow(sample.data)) 
        else group <- sample.data$group
      if (is.null(names(h))) names(h) <- sort(unique(group))

    regular.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known"))
    augment.fit <- try(rrreg(y ~ x1 + x2, data = sample.data, p = 0.5, p0 = 0, p1 = 0.5, design = "forced-known", 
      h = h, group = group, matrixMethod = weighting))
    }
    
    sample.X <- cbind(1, sample.data$x1, sample.data$x2) 
    sample.0 <- sample.X[which(sample.data$x1 == 0), ]
    sample.1 <- sample.X[which(sample.data$x1 == 1), ]

    result <- rbind(
                  cbind(mean(logistic(sample.0 %*% coef(regular.fit))), 
                      mean(logistic(sample.1 %*% coef(regular.fit))), 
                      0), 
                  cbind(mean(logistic(sample.0 %*% coef(augment.fit))), 
                      mean(logistic(sample.1 %*% coef(augment.fit))), 
                      1)
                  )
    colnames(result) <- c("x1.0", "x1.1", "augment")

    result

 }  

